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Abstract 

In a two-part publication, we propose and analyze a formal foundation for practical calcula- 
tions of vibrational mode lifetimes in solids. The approach is based on a recursion method 
analysis of the Liouvillian. In the first part, we derived the lifetime of vibrational modes in 
terms of moments of the power spectrum of the Liouvillian as projected onto the relevant 
subspace of phase space. In practical terms, the moments are evaluated as ensemble aver- 
ages of well-defined operators, meaning that the entire calculation is to be done with Monte 
Carlo. In this second part, we present a numerical analysis of a simple anharmonic model 
of lattice vibrations which exhibits two regimes of behavior, at low temperature and at high 
temperature. Our results show that, for this simple model, the mode lifetime as a function 
of temperature and wavevector can be simply approximated as a function of the shift in fre- 
quency from the harmonic limit. We next compare these calculations, obtained using both 
Monte Carlo and computationally intensive molecular dynamics, with those using the lowest 
order moment formalism from the Part I. We show that, in the high-temperature regime, 
the lowest order approximation gives a reliable approximation to the calculated lifetimes. 
The results also show that extension to at least fourth moment is required to obtain reliable 
results over a full range of temperatures. 
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1. Introduction 

The calculation of vibrational mode lifetimes using the Green-Kubo relation pQ in solids 
is, at present, a computationally expensive process. The formalism presented in Part I of 
this work [2] provides an approach to numeric computation of these lifetimes utilizing only 
ensemble statistics, allowing for much quicker computation. The strength of this formalism 
is greatly enhanced as it does not require the use of computationally expensive molecular 
dynamics needed to calculate autocorrelations of occupation numbers in the solid up to large 
times. Rather, mode lifetimes can be approximated by examining the properties of a given 
system over the ensemble, which requires no molecular dynamics (MD). 

The mode lifetime is defined here by the Green-Kubo relation 

f+oo 

r k = dt X k(t) (1) 



where r k is the lifetime of mode k, and Xkif) is the mode autocorrelation function 

(6n k (0) Sn k (t)) 

Xkit) = wjm (2) 

Here, Sn k is the fluctuation of the occupation number of a given mode, and the angular 
brackets indicate an average over the equilibrium distribution. The straightforward but te- 
dious method of calculating r k requires evolving a given state forward in time using MD, 
calculating x(t) and averaging over the equilibrium ensemble. The improved method, in- 
troduced in Part 1, is based on applying the recursion method [3j HJ [5] to the Liouvillian 
operator [HI Ej, L, defined by 

f dH d dH d . 

with q and p the coordinate and momentum variables of the system and H the Hamiltonian. 
By successive applications of the Liouvillian on 5n k we can generate a sequence of orthonor- 
mal functionsjJJ |8]. With this sequence as a basis, the Liouvillian takes a special tridiagonal 
form. Using this scheme, Xkif) can be simply related to the resolvent R(u) of the Liouvillian 
(defined by R(u>) = (u> — L)^ 1 ) . The resolvent can then be expressed as a continued fraction. 
More precisely, the auto-correlation is related to the projection of the resolvent onto 5n k . It 
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can then be shown that can be expressed in terms of the moments of x(w)> t ne Fourier 
transform of x{t)- The power of the method comes if it is possible to obtain a reasonable 
approximation from only a few low-order moments. 

In order to evaluate the effectiveness of this method, we compare the lowest-order calcu- 
lation to numerically exact results for a simple model of anharmonic lattice dynamics. To 
this end, we have determined tj, two ways, first, by calculating it directly using MD and 
Monte Carlo, and second, using the recursion method described above, truncated to lowest 
moment. We also characterize some of the interesting behavior of the model used. 



2. Methods 

Our model is based a continuous vector-like quantity defined on a three-dimensional 
lattice. The underlying lattice structure is simple cubic (8x8x8, with 512 sites) with periodic 
boundaries. Nearest neighbor lattice sites are connected by anharmonic potentials such that 
the Hamiltonian of the system is given by 

H = J2l\&\ 2 + E Viji-dj) (4) 

i <i,j> 

where pi is the momentum of the particle at lattice site, i, and V(d) is given by 

v{d) = \\A 2 + ^ (5) 

The coupling occurs for every pair which are nearest-neighbors. 

Note that this model describes a vector degree of freedom (a) on each lattice site, which 
is connected to neighboring degrees of freedom. However, it does not capture particle mo- 
tion. In other words, since the neighbors in the lattice are fixed, there is no mass flow, but 
only vibrations about a fixed lattice. In this model, there is no melting transition, which we 
consider to be an advantage because we can examine the range of validity of our calculations 
over a virtually unlimited range of temperatures. This is in contrast to the work of Ladd, 
et al. [9] who studied mode lifetimes in a Lennard- Jones system, which does have a melting 
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transition, so that their studies were effectively limited to a narrower range of temperatures. 
Furthermore, by choosing a Hamiltonian which does not include particle flow, we are able 
to selectively study the nature of anharmonic lattice vibrations. 

We also note that this model displays different behavior in two temperature regimes. At 
low temperature, the system behavior is dominated by the harmonic part of the potential, 
with only a weak anharmonicity evident in the dynamics. This is characterized by a har- 
monic heat capacity and long lifetimes for the modes. At higher temperature, the system is 
dominated by the quartic part of the potential, which is characterized by a different heat ca- 
pacity and shorter lifetimes for the modes. The transition between the two regimes appears 
to be at a temperature of about 10. 

In order to calculate Xk{t), an initial position is chosen in phase space, using a Monte 
Carlo sampling at some fixed temperture, and this state is propagated forward in time us- 
ing the velocity Verlet algorithm. At regular time intervals, x(t) can be calculated from 
the current and initial states, following the canonical transformation laid out in Part I to 
determine the occupation number, n. As is noted there, we can minimize the fluctuation in 
x(t) by carefully choosing the frequency used in the transformation, effectively choosing the 
anharmonic frequency, u a as our transformation frequency. This can be done by choosing 
the frequency which maximizes r^, as is explained in Part I. x(t) i s then calculated as an 
average over fifty-thousand initial Monte Carlo points. 

In part I, we showed how r(k) could be related to the moments (/i n ) of the resolvent, 
suitably projected onto the appropriate function in phase-space corresponding to the mode 
with wavevector k. That is, we obtained the exact result: 



r(k) = F(/i 2 ,/i 4 ,/i 6 ,...) 



where the moment can be calculated directly from the Liouvillian by 




((^ fc ) 2 ) 
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This can be re-expressed using dimensional analysis [TU] in the form 

r{k) =t 2 (/c)F( T4 , 76 ,...) 

where 




and the 7 n are dimensionless forms of higher moments, such as 74 = fii/i^) 2 - Roughly 
speaking, the second moment describes the width of the appropriate part of the power spec- 
trum is the Fourier Transform of x{t))i while the 7's describe the shape of that power 
spectrum. 

The approximation that naturally suggests itself here is to assume that whatever changes 
in the power spectrum occur with temperature, they can be tracked by just a small number 
of the lowest moments of the power spectrum (or equivalently, the second moment and a 
small number of 7's). The simplest form of this approxmation would be to rely entirely on 
the second moment (^2); which is equivalent to assuming that the 7„ are all independent 
of temperature (that is, the power spectrum only changes width but not overall shape). 
There has been considerable success in other types of physical systems with this moment- 
approximation, most notably in the calculation of electronic structure of materials [8| [TT| [12] . 

We therefore evaluate here the second moment approximation, for which r(k) oc T2(k) 
where the proportionality constant is some function of the 7n 's, which are all assumed to be 
independent of temperature. We apply two tests. First, we calculate r(k) according to the 
Green-Kubo prescription (which involves MD), and the T 2 (k) using Monte Carlo alone. If 
the second moment approximation is successful, the ratio ^^(/c) should be independent of 
wavevector and temperature. The second test is more stringent. If the time-scale of mode 
decay is simply proportional to r 2 , then this should be reflected in the auto-correlations 
x{t)- Namely, if we scale the time variable in auto-correlation function for each mode by 
r 2 , we should see a "data collapse". In other words, plotting x(t/ T 2) for all wavevectors and 
temperature should yield a universal curve. The detailed shape of that universal, scaled auto- 
correlation is determined by the shape parameters (the 7 „'s). The "data collapse" test is 
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Figure 1: A typical autocorrelation, for our model 



therefore a more rigorous determination of the validity of the second-moment approximation. 



As we will show in the following, the high-temperature behavior is quite well-represented 
by the second moment approximation, but that the lower-temperature behavior shows signif- 
icant deviation. This indicates that more moments may be needed to accurately approximate 
the low temperature behavior. 



3. Numerical Results 

In this section, we review first the results of the numeric calculation of the Green-Kubo 
lifetimes for our model Hamiltonian. A sample of the autocorrelation as a function of time is 
given below (Fig. fl]). It is first interesting to note that, due to the finite size of the ensemble 
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Figure 2: r as a function of temperature in the [220] direction. The solid line represents r oc T 1 ^ 3 . Note 
that for temperatures less than 10, the slope is steeper. 

used, the autocorrelation does not go strictly to zero for large times, so the termination 
of the autocorrelation is nontrivial. It would appear that, because of the finite number of 
degrees of freedom, the ensemble retains some memory of its initial state for arbitrarily large 
times. An examination of these residual "tails" as a function of system size shows that as 
the system becomes large, the magnitude of the tails quickly goes to zero. For our finite 
systems, the r integral was truncated after the autocorrelation dropped below 0.5 percent 
of the original value. 

We first show how r varies with temperature for a few wavevectors and note that the 
general behavior is the same for all of them (Fig. [2]). As can be seen, the dependence on 
temperature to separate into two regions, high and low temperature being above and below 
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about 10. For high temperatures, r is related to temperature by a simple power law, namely 

r(T) oc (6) 

This is different from the behavior discussed by Ladd for Lennard-Jonesium, namely 
r(T) oc T _1 . However, they expect this behavior only at low temperature. We have not 
determined the behavior of our model at low temperatures, because the computational times 
increase substantially owing to the extended lifetimes at low temperature. For the remain- 
der of this paper, we will focus on the high temperture regime and will only note apparent 
deviations as we approach the transition to low temperature. 

Next we examine r along the (100), (110), and (111) directions (Fig. [3]). We see here 
that, as expected, r decreases with increased wavevector (as well as increased frequency), 
but because of the small size of the system, little more can be observed. 



In order to see more directly how r is determined by the temperature and wavevector for 
this Hamiltonian, it is instructive to examine the behavior of the anharmonic mode frequency 
(u a (k)). For our model Hamiltonian, we find a simple relation between the anharmonic fre- 
quency, the temperature, and the frequency in the harmonic limit (uo(k)). The harmonic 
frequency is of the mode with wavevector k is 



uj (k) 



2nk x 2irk v 2irk z 
,6 — 2 cos — 2 cos — 2 cos — - — (7) 

\ <>x <"y I'Z 



As can be seen in Fig. [4J the ratio between the frequency shift (5uj(k) = oj a (k) — uJo(k)) and 
Co>o, depends only on temperature (that is, not on wavevector) , and is well approximated by 
the power law 

^|=/(T) = (T/T„)l (8) 

While the exponent and possibly the simple form of the relationship is specific to our Hamil- 
tonian, it is still an intriguing result. Based on the similarities in the exponent in Eq. [6] and 
[8j one reasonably suspects a relationship between the 8co and r. This suspicion is confirmed 
by examining r as a function of 5co (see Fig. |5l) for all wavevectors at a constant temperature, 
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Figure 3: Tau as a function of wavevector in the (a) (1,0,0), (b) (1,1,0), and (c), (1,1,1) directions 
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Figure 4: Suj/uj as a function of temperature. This ratio is the same regardless of the wavevector examined. 
The solid line represents Suj/uj oc T 1 / 3 

where we see again a power-law dependence. 

The implication of this relationship is that, at a fixed temperature the lifetime of a given 
mode depends dominantly on the anharmonic frequency. This frequency, in turn, depends 
on two independent factors, the harmonic frequency and the temperature. 

Knowing the mode lifetimes as a function of wavevector and temperature through our 
numerical simulation, we can now turn to the approximate methods discussed in the com- 
panion work. 
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Figure 5: r as a function of Sui at a temperature of 2000 for all nondegenerate wavevectors of the model 
used. 
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4. Approximate Methods 

Because its determination does not require the use of any molecular dynamics, r 2 can be 
calculated with much less computational expense then the numerical approach considered 
above. Generally we have found that 200,000 MC steps is sufficient to converge r 2 for all 
temperatures and wavevectors considered here. 

If we first examine the ratio 7&/T2 as a function of temperature for a particular wavevec- 
tor, we see that, in the high temperature regime, the ratio quickly converges to constant 
value (Fig. [7]). While only one wavevector is shown, note that this behavior holds for all 
wavevectors. This is expected, due to the separable nature of the anharmonic frequency 
for our Hamiltonian and the dependance of r on this frequency mentioned above. Because 
t can be determined from the shift in frequency, the rest of them should display the same 
behavior. Because the individual wavevectors may have different shapes, the ratio may con- 
verge to different values. However, it remains true that this value is constant for a given 
wavevector for sufficiently high temperatures. The implication is that only r 2 is changing 
with temperature, and the 7 factors discussed above are constant. 

With this notion that the shape (as opposed to the width) of x(t) is invariant, we also 
test whether rescaling the elapsed time with respect to r 2 will cause the data to collapse 
onto the same curve. This is a more stringent test than the previous one as it indicates that 
x{t) and not just r can be rescaled with r 2 . Note, in fact, that the collapse of x{t) implies 
the convergence of Tfc/r 2 seen in Fig. [6] As Fig. [7] shows, this is clearly the case in the high 
temperature regime. In this limit, approximating r at the level of the second moment is 
sufficient to determine its behavior as a function of temperature. 

Because this agreement begins to break down at lower temperatures, it may be neccesary 
to include higher moments as the shape of x(t) changes from the high temperature limit. As 
the shape of the autocorrelation changes, some number of the 7 factors are also changing. 
By calculating them, and taking into account their effect on the system, a more accurate 
prediction of r should be possible. 
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Figure 6: t/t2 as a function of temperature. We see that ti is sufficient determine r for high enough 
temperatures. 



5. Conclusions 

By analyzing the mode lifetimes, through in an anharmonic lattice, with a sim- 

ple, quartic interaction between nearest neighbors, in three dimensions, we can draw several 
important conclusions. First, for this particular Hamiltonian, the lifetime of a given mode 
depends in a simple way on the anharmonic frequency of the mode, which in turn depends 
separably on the harmonic frequency and the temperature. Both of these relationships would 
seem to be peculiar to our model Hamiltonian, but they allow a simplified analaysis of the 
results. Second, seen as a function of temperature and frequency directly, r shows a simple 
behavior. More generally, and most central to this work, the methodology developed in Part 
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Figure 7: The data collapse of x(t) for the [222] wavevector at several different temperatures. Because the 
data approach a universal curve, %(t) and therefore r are well-approximated, up to a constant, only by T2- 
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I is shown to be a reliable convergent method of determining the mode lifetimes with much 
lower computational cost than any present method of which the authors are aware. While 
deviations are observed for low temperatures, related to the form of the potential used, using 
only the first non-zero moment, x(t) collapses effectively onto a universal curve, and r is 
correctly approximated at high temperatures up to a constant to within a few percent. 

While the Hamiltonian used here is a relatively simple one, which might be analyzed ef- 
fectively using more direct methods, the scope of this new method should not be understated. 
Since calculation of lifetimes using moments requires no molecular dynamics, only averages 
over the ensemble, the computer time required for a calcuation is reduced dramatically, in 
the authors' experience, by at least an order of magnitude or more at high temperatures, 
where jtau is small and even more for lower temperatures as jtau increases and molecular 
dynamics calculations must go to higher times. Furthermore, the method can be readily 
extended to more complicated systems. As long as 5n can be found for a member of the 
ensemble, r 2 can be derived. In addition, this formalism provides a new way of examining 
vibrational mode lifetimes which may provide new insights into their behavior. This im- 
proved method offers not only significant advantages in both improving the methodology of 
calculation but also possible insights into the mechanism of dissipation. 
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